import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import make_blobs
import time

def kmeans(dataset,cluster_n,):
    #任取3个元素作为其的中心选取cluster_n个，索引在len(dataset)内

    #创建全局最优解变量 
    best_center=None
    best_label=None
    best_avdistance=np.inf
    best_distance=None

    for i in range(10):
        # 3*2 (cluster_num*xy)初始中心点
        initial_center=dataset[np.random.choice(len(dataset),cluster_n,replace=False)]
        
        center=np.zeros((cluster_n,2))

        #  100*3 (sample*cluster_num) 每个点距离三个中心点的距离
        distances=np.zeros((len(dataset),cluster_n))

        #只要原始数组和变换后的数组不同就不断更新中心值
        while(np.sum((initial_center-center)**2)>1e-8):

            #加入判断，若不为0矩阵才进行赋值操作。
            if not(center==0).all():
                initial_center=np.copy(center)

            #计算距离
            for j in range(cluster_n):
                distances[:,j]=np.sqrt(np.sum((dataset-initial_center[j])**2,axis=1))#sum将x**2和y**2相加
        
            #获取最小值的索引
            # 100*1 (sasmple*choice)每个点属于的标签值
            label=np.argmin(distances,axis=1) 
    
            for j in range(cluster_n):
                #initial_center[j] = np.mean(dataset[label == j], axis=0)也可已使用mean直接求
                #必须限定方向，否则会自动算这个矩阵的平均值
                center[j]=np.sum(dataset[label==j],axis=0)/np.sum(label[:]==j)

        #获取所有点距离自己最近点的平均距离
        for j in range(cluster_n):
            distances[:,j]=np.sqrt(np.sum((dataset-center[j])**2,axis=1))

        #计算每个点距离自己最近的簇的距离的平均值  
        #距离最近簇的距离
        distance=np.min(distances,axis=1)

        #平均距离
        average_dis=np.mean(distance)
        
        #迭代寻找最优簇的点
        if(average_dis<best_avdistance):
            best_avdistance=average_dis
            best_label=label
            best_center=center
            best_distance=distance
        #手动设置异常点的判断条件(distance>2*average_dis)
        #创建了一个bool数组使用bool索引来创建异常点
    outliers=dataset[best_distance>2*best_avdistance]
    print("异常点为：",outliers)
    
    '''for i in range(100):
        #超出两倍即为异常点
        if distance[i]>2*average_dis:
            print(dataset[i])
    '''
    
    # 可视化
    plt.scatter(dataset[:, 0], dataset[:, 1], c=best_label, cmap='viridis')
    plt.scatter(best_center[:, 0], best_center[:, 1], c='red', marker='x')
    plt.scatter(outliers[:, 0], outliers[:, 1], c='black', marker='o', s=100, alpha=0.5, label='Outliers')
    plt.legend()
    plt.show()

# 生成数据
X, y = make_blobs(n_samples=100, n_features=2, centers=3, random_state=int(time.time()))

# 将数据点和标签转换为NumPy数组
dataset = np.array(X)#  100*2(sample*xy)
kmeans(dataset,10)
